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DIRECT-NUMERICAL AND LARGE-EDDY SIMULATIONS OF A NON-EQUILIBRIUM 

TURBULENT KOLMOGOROV FLOW* 

S.L. WOODRUFF* , J.V. SIIEBALIN* , AND M.Y. HUSSAINI$ 


Abstract. A non-equilibrium form of turbulent Kolmogorov flow is set up by making an instantaneous 
change in the amplitude of the spatially-periodic forcing. It is found that the response of the flow to 
this instantaneous change becomes more dramatic as the wavenumber of the forcing is increased, and, at the 
same time, that the faithfulness with which the large-eddy-simulation results agree with the direct-numerical 
results decreases. 

Key words, turbulence, modeling, large-eddy simulation 

Subject classification. Fluid Mechanics 

1. Introduction. How, and by how much, one must modify present steady-state sub-grid models in 
order to perform largo-eddy simulations of non-equilibrium turbulent flows is very much ail open question. 
The ever-growing computational resources at our disposal provide the opportunity for applying large-eddy 
simulations to non-equilibrium as well as other increasingly-complex flows, but there is no guarantee that the 
sub-grid models currently in use will be up to the task of adequately represent ing the effect of the sub-grid 
scales on the resolved flow. In the case of the non-equilibrium flows of interest hen 1 , for which the statistical 
state of the turbulence varies with time, there is clearly a failure of a steady-state mock 1 ! one which bases 
its prediction of the Reynolds stress on the current, instantaneous, resolved velocity field to reflect the 
finite time lag inherent in the reaction of the actual sub-grid dynamics to temporal changes in the resolved 
field. The extent to which this effect actually affects the large-eddy simulation results is another question, 
and clearly one strongly influenced by the rapidity of the resolved-field variations. The present investigation 
seeks to answer this question for one simple flow. 

The choice of a turbulent flow to use in the study of this question is not a trivial one, given our desire 
to work with a flow whose computation does not. require too much computer time, whose numerical analysis 
does not require especially-sophisticated numerical techniques and whose physics bears the closeset possible 
relationship to turbulent, flows of real technological interest. Homogeneous, isotropic, turbulence is the 
simplest choice, but the applicability of results for this flow to real flows is problematic, at best. Homogeneous 
shear flow is a possibility, but the additional complexity of the numerical techniques required to deal with 
the cumulative shearing regridding, [8, 2, 9, 11]) make it an unappealing one. Simulations of even the 
simplest, of those turbulent flows realizable in the laboratory, such as channel flow, are both computationally 
intensive and require special attention near walls, both to the numerics and to the modelling. 

Shebalin and Woodruff [14] have 1 advocated the study of turbulent. Kolmogorov flow as a way to learn 
about, turbulent shear flows without the excessive numerical difficulties and computational expense of the 
flows referred to above. This flow, in which fluid in an infinite domain is driven by a periodic body force. 

* 1 his research was supported by the National Aeronautics and Space Administration under NASA Contract No. NASt- 
19180 while the first and third authors were in residence at the Institute for Computer Applications in Science and Engineering 
(ICA8E), NASA Langley Research (’enter, Hampton, YA 2:1081-2199. 

* School of Computational Science and Informat ion Technology. Florida State University. Tallahassee, FI. 32300-1120 

* Aerodynamic and Acoustic Methods Branch. NASA Langley Research Center, Hampton. YA 23081-2199; currently at 
Johnson Space Flight Center 

& School of Computational Science and Information Technology, Florida State University, Tallahassee. FL 32300-1120 
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was proposed by Kolmogorov as a model problem for the study of stability issues, and a significant body 
of literature addresses the stability aspects of this flow (see, for example, [16], and the references therein). 
Additionally, there has been some work in the numerical simulation of turbulent Kolmogorov flow in two 
dimensions [13, 12] and in three dimensions [14]. 

In addition to the computational and numerical advantages afforded by the Kolmogorov flow, it has 
a further advantage for the testing of turbulence models in that there is the potential for allowing many 
different turbulent flow features to be studied, due to the freedom to choose the body force as one pleases. 
In the present work, use is made of this freedom to investigate the response of the flow to instantaneous 
changes in the amplitude of the forcing at different shear magnitudes. 

A number of approaches have been proposed for the development of time-dependent models for non- 
equilibrium flows. Yakhot and Smith [15] presented a time-dependent eddy viscosity for k-t models; a more 
general expression of the idea behind this model provided the basis for the stress-relaxation model of Yakhot 
et al [18]. The two-scale Direct Interaction Approximation approach to the development of turbulence 
models proposed by Yoshizawa [20] was employed by Yoshizawa and Nisizima [21] to derive a time-dependent 
model; Rubinstein [10] also used Kraichnan’s DIA [4], to motivate time-evolution equations for the Reynolds 
stresses. The general approach to developing turbulence models examined in [17], based on the fundamental 
idea behind Yoshizawa’s work, was used to derive a history-integral model for non-equilibrium turbulent 
flows, independently of any time-scale-separation assumptions. Two simplified versions of this model were 
tested numerically in large-eddy simulations of plane-channel flow by Nwafor and Woodruff [6]; they found 
that the large-eddy simulations with the Smagorinsky model did in fact fail to reproduce the DNS results, 
and the new models did improve the results in several respects. 

Reported in this paper are results of a comparison test for Kolmogorov flow similar to that of Nwafor 
and Woodruff. Simulations for Kolmogorov flows with spatially-sinusoidal forcing at three different wave 
numbers are set up and run until a steady state is reached; then the forcing amplitude is doubled and the 
relaxation of the flow to a new steady state is examined. These three non-equilibrium turbulent Kolmogorov 
flows were solved by direct numerical simulation and large-eddv simulation with the plane-averaged history- 
integral model (a simplified model proposed in [6]), with the dynamic sug-grid scale model [7] and with the 
Smagorinsky model using two values of the Smagorinsky constant: one consistent with the predictions of the 
history-integral model (C*) and the other consistent with the predictions of the dynamic model (C f /). (The 
latter is different for the three values of forcing wave number, the former is not.) 

In contrast to the findings of [6], it was found that there was little difference between the results of the 
LES with the new model and those of the LES with the Smagorinsky model employing CV What differences 
there were increased as the 1 wavenumber of the forcing increased. The LES reproduced the DNS results fairly 
well at the lowest wave-number forcing, but became progressively worse as the wave number of the forcing 
increased. The dynamic model was an improvement over the Smagorinsky model with C/, and the history- 
integral model at the higher- wave-number forcings, but it, too, exhibited discrepancies when compared with 
the DNS. The Smagorinsky model with C t i did nearly as well as the dynamic model, indicating that the 
primary advantage of the dynamic model in this application is its automatic initial determination of the 
Smagorinsky constant; once the constant has been determined near the beginning of the simulation, there 
seems to be little advantage gained from the dynamical model's ability to further modify the constant in 
time or space. 

The specifics of the physical problem to be solved are given in the following section. Section 3 gives 
information about the models employed in the large-eddy simulations and Section 4 is concerned with 
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numerical details, m particular, details of the implementation of the new model. The results are described 

in Section 5 and the implications of these results for the modelling of non-equilibrium turbulent flows are 
discussed in Section 6. 

2. Description of the Flow and Its Properties. The Kolmogorov flow to be considered in this work 
is contained within a periodic box whose sides are of length 2tt. Energy is input by an external, artificial, 
deterministic, body force which is introduced into the equations on the right-hand side. We use a Cartesian 
coordinate system, with axes x, y and c. 

The external forces to be used throughout this work involve a single non-zero component in the .r direction 
which depends only on the 5 coordinate. As a result, the flows to be examined are parallel, as are Ornette 
and Poiseuille flow, m the sense that the mean motion is in one direction and varies in a perpendicular 
direction. The flow corresponds to Kolmogorov flow if the non-zero component of the force is sinusoidal in c 
and we shall here consider only sinusoidal forcing at various wavenumbers. We are thus interested in solving 
the Navier-Stokes equations with a force of the form f = kjv^mkfz t inserted on the right-hand side. The 

chaiactciistk velocity r 0 and the forcing wavenumber k/ may be used to nondimensionalize the equations 
of motion: 


(2.1) 


DU V7 , 1 T-> 

~Dt ~ ~ Vp + u + sni ; 


the Reynolds number Re is Vo/kfis. 

In addition to the parallel nature of the flow, we may conclude from the fact that nothing depends on 
the y coordinate that the flow properties must not vary under an inversion of the y coordinate (y -y- 
the origin is at the midplane of the box.) This implies, for example, that all Reynolds averages involving 
the y component, of velocity must be zero. The only nonzero off-diagonal term of the Reynolds stress tensor 

is thus the < uw > term, and consequently this term is the one which produces energy for the turbulent 
velocity field. 


3. Description of the Models. In this section, we briefly review the reasoning behind the Smagorin- 
sky model, examine the issue of whether a time-dependent model is necessary for large-eddy simulations of 
non-equilibrium turbulent flows, and review the derivations of the history-integral model [5] and the plane- 
averaged history-integral model [6], a simplification which allows us to incorporate the history integral into 
a calculation in a computationally-feasible manner. Finally, we describe the application of the dynamic 
sub-grid scale model to this flow. 

The well-known Smagorinsky model may be motivated by an appeal to the isotropy of the modeled 
subgrid scales, which leads to the form 

R’j = + l 'rS,j\ 

k being the kinetic energy of the subgrid scales, the turbulent eddy viscosity and S u the rate-of-strain 
tensor of the resolved velocity field f/,: 


(3.2) 


Sn = 




One then has to determine the form of the eddy viscosity, v r , and dimensional analysis tells us that once 
we have chosen to assume that the modeled subgrid scales depend only on the dissipation rate (following 
Kolmogorov s analysis [3]) the eddy viscosity must have the form u, = const. r-V». The dissipation rate for 
the subgrid scales is determined by assuming that the energy balance for the subgrid scales is dominated by 



the local production RijSij and the local dissipation e. If this assumption is valid, then the kinetic energy 
equation for the subgrid scales reduces to a statement of the equality of the local production and the local 
dissipation. This relation may be used along with the Smagorinsky model itself (3.1) to eliminate e, giving 
an expression in terms of the rate-of-strain tensor only: 

(3.3) Rjj = - j + (C s A)-(S m nSm») ^ S>j 

Here A is a length scale characterizing the grid size employed in the calculation and the constant is determined 
empirically by fitting with experiments or direct-numerical simulation data; its value is typically taken to be 

0 . 1 . 

The crucial point about the Smagorinsky model for the present discussion is that it is based on scale- 
separation assumptions, both in space and in time. Gradient-transport models in general result from scale- 
separation assumptions: for example, the Newtonian-fluid viscous-stress relation is derived in the Chapman- 
Enskog formulation by taking advantage of the large separation in scales between the continuous fluid motion 
and the discrete molecular motion. Formal derivations of the Smagorinsky model also employ such a scale 
assumption (see, for example, [20, 17]). 

Since there is in turbulence no spectral gap — a range of scales between the largest and the smallest in 
which the energy is neglegible the scale-separation assumption is fundamentally not valid and we have to 
ask whether it should be removed. Certainly the success of the Smagorinsky model in simulations of many 
different types of flows suggests that it at least, provides a useful formula, independently of the invalidity 
of the scale-separation assumption. The question of concern is really whether or not there are flows which 
are simulated better by a model which does not embody this scale separation assumption then they are 
simulated by the Smagorinsky model. 

A model which removes the temporal scale-separation assumption has been derived by Woodruff [17], 
and has been applied to the non-equilibrium turbulent flow problem of accelerated plane-channel flow by 
Nwafor and Woodruff [6]. There it was found that there are in fact significant features of that flow which 
compare better with direct numerical simulatioins when simulated with the new model than when simulated 
with the Smagorinsky mdoel. Thus, there is some evidence that, from a purely practical standpoint, there 
are situations for which a time-dependent model does improve the results of large-eddy simulations of non- 
equilibrium turbulent flows. In the present paper, we examine whether this is true for Kolmogorov flow. 

The time-dependent model used in the present investigation was derived in [17], by employing a fairly 
general technique for the derivation of turbulence models using ideas from the analytical theory of turbulence. 
The fundamental idea behind the approach goes back to Yoshizawa [20] and even further to Grow [1] and is a 
sort of rheological approach to turbulence wherein one separates the motion into the modeled portion and the 
resolved portion and treats those portions as distinct entities. A solution is derived for the modeled portion 
which involves functions characterizing the resolved portion; this solution may be used to give expressions for 
the modeled contributions to the resolved portions of the flow. Specifically, for the modelling of turbulence 
we divide the flow into the resolved portion and the subgrid, modeled, portion and derive a perturbative 
solution for the subgrid velocities, substitute that perturbative solution into the definition of the Reynolds 
stress and so derive an expression for the Reynolds stress which may be inserted into the Reynolds-averaged 
equations (or, in the case of large-eddy simulations, the filtered equations). The particular approximate 
solution for the subgrid velocities used in [17] was based on the assumption that if the commonly held 
assumptions about the universality of turbulence at the small scales are valid, we can approximate the sub- 
grid scale velocities in terms of a dominant part which is the contribution to the velocity from the universal 
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component and depends (according to the assumption of universality) only on a small number of descriptors 
of the resolved flow, such as the dissipation rate. Corrections to this dominant term arise from more specific 
details of the resolved flow, such as the rate-of-strain tensor. Thus, we have what is essentially an expansion 
in the magnitude of the rate of strain. The terms in this expansion may in principle be computed by any 
means that is available and convenient, using the theory of turbulence, or even using numerical or laboratory 
experiments. The models developed in the examples of Woodruff [17] were derived on the basis of turbulence 
approximations derivable from the Direct Interaction Aproximation of Kraiclman [4], The one employed in 
the present numerical calculations, as well as the calculations of Nwafor and Woodruff [6], is based on the 
e-expansion renormalization-group results of Yakhot and Orszag [19], In fact, any specific theory used to 
make a computation along these lines would lead to a history-integral model of the type discussed below, 
the only differences would be in the details of the kernal function and it is likely that these differences would 

not affect the results of a tubulence model calculation drastically. The reader is referred to [17] for details 
of this derivation. 

It, is helpful ill understanding the assumptions underlying this new model to sketch a more intuitive 
derivation. In this derivation, we make several assumptions and then ask what is the most general form an 
expression for the Reynolds stress may take consistent with those assumptions. The assumptions are that the 
Reynolds stress is linear in the strain rate, that the model is isotropic and and that the model is, spatially, a 
gradient transport model. The consequences of these three assumptions are that the proposed model must 
be identical in form to the Smagorinsky model, as far as the spatial dependence goes. The difference lies in 
the temporal dependence, which may be represented without loss of generality as a linear integral operator: 

r i 


(3.4) 




2, . f 1 

- 3 ^ + l 


A'(M) Sjj(s) (Is: 


for clarity, we have shown explicitly only the time dependence of R u and S ir 

The kernal A'(t, s) is determined by the analysis of [17] in terms of the response function and correlation 
function of the modeled subgrid velocity held and the expressions derived there will be used in the calculations 
to follow. It is, however, possible to argue from dimensional considerations that the kernal has the form 


(3.5) 


A '(/,*) = 




where k, is the wavenumber of the smallest resolved eddies in the numerical calculation and t, is the 
corresponding eddy turn-over time. It is natural to expect that, more recent history is more important than 
less recent history, so the non-dimensional function F should be a monotonically decreasing function of its 
argument. The analysis of [17] gives 


(3.G) 


F(.r) 


= '[ 


~ r 


-2x 


+ 2xEi(-2x)< 


and this expression will be used in the present calculations. 

If the history-integral model were simplified on the basis of a scale-separation assumption in time, it 
would reduce to the Smagorinsky model as a sort of steady-state limit. (That is, the rate-of-strain tensor 
is assumed to be essentially constant as far as the evaluation of the history-integral is concerned.) This 
calculation was performed in [17], where the predicted value of the Smagorinsky constant was found to be 
quite close to the value traditionally used in channel flows. 

The computations required to actually compute this history-integral expression for the Reynolds stress 
m the course of a numerical calculation would lie quite extensive, and certainly prohibitive for a practical 



calculation. For this reason, we choose to take advantage of the homogeneity of the Kolmogorov flow under 
consideration here in the x and y directions to consider the possibility that most of the history effects may 
be captured by considering only the time history of the rate-of-strain tensor averaged over * - y planes. 
This approach was also employed in the plane-channel calculations of Nwafor and Woodruff [6]. Such a 
plane-averaged approach may not be implemented rigorously, because it is not possible to rewrite equation 
(3.4) in terms of only plane-averaged quantities without breaking up some averages of products of quantities 
into products of averages. Once this is agreed to, however, and once the basic Smagonnsky model is added 
and subtracted from the history-integral model, the plane-averaged form of the history-integral model may 
be written 

* dt' K p (t - t') < > P - < v{t) >,,< Sij > p 

(The notation < • > p indicates plane averaging and the subscript on the kernal indicates that the kernal 
is computed using the plane-averaged dissipation rate.) This expression has been specifically constructed 
so as to be the basic Smagorinsky model plus a history-integral correction term, to aid in the numerical 
implementation of the model. The implementation will be sketched in the following section. 

The dynamic sub-grid scale model is implemented following [7]. A test filter is defined which filters the 
resolved velocity field such that the smallest resolved length and time scales are eliminated. One then seeks a 
value for the Smagorinsky constant which yields the same values for the modeled Reynolds stress from both 
the resolved velocity field and the test-filtered velocity field. Given that there are six independent components 
of the stress tensor and but one Smagorinsky constant , this requirement leads to an overdetermined system 
of equations; of the several approaches to finding the “best” value for the Smagorinsky constant under these 
constraints [7], that proposed by Lilly [5], where the square of the residual error is minimized, is employed 
here. Additionally, as is common in this sort of parallel-flow shear-layer problem, the computation of the 
Smagorinsky model is performed with quantities averaged over x-y planes in order to ensure numerical 

stability. 

4. Numerical Considerations. The DNS code used in the present work is that, of Shebalm and 
Woodruff [14]. The code is spectral with Fourier modes in all three spatial directions and incorporates 
the viscous term implicitly in time-stepping based on a predictor-corrector algorithm. The nonlinear terms 
are handled by fast-Fourier transforming back into physical space, multiplying and then transforming into 
Fourier space. It was found that a resolution of 64 3 was sufficient for an accurately-resolved solution for the 
parameters employed in this study. 

In order to perform large-eddy simulations of the Kolmogorov flow, subroutines were added to the DNS 
code which compute the Reynolds stress based on the Smagorinsky model, the dynamic model and the 
plane-averaged history integral model. For the runs described here, the original time-stepping of the DNS 
code was retained, where only the (molecular) viscous terms are treated implicitly: the entire Reynolds stress 
term is treated explicitly. Some experimentation was done with a time-stepping algorithm incorporating the 
eddy viscosity averaged over x - y planes into the implicit part of the computation, but the second-order 
differences employed in the 2 direction in this altered algorithm were found to degrade the accuracy of the 
calculation over long times at higher values of shear and this algorithm was abandoned. 

The Reynolds stress for the Smagorinsky model is computed by calculating the rate-of-strain tensor for 
the resolved velocity field and then computing the Smagorinsky approximation for the stress according to 
(3.3). The dynamic model is computed similarly, with the addition of the computation of the Smagorinsky 

constant itself. 


(3.7) 


R,j — 


u(t)Sij + 


Implementation of the plane-averaged history integral model is of course more complex, since it is 
necessary to store the history of the rate-of-strain tensor and compute its integral. As described in Section 
3, we implement the plane-averaged history integral model by inserting it as a correction factor to the 
Smagonnsky model so the routines described above for the Smagorinskv model remain intact. In addition, 
we add routines to calculate the plane-averaged rate of strain tensor, the kernal of the history integral and the 
history integral itself. The dissipation for use in the computation of the history integral kernal is computed 
according to the local production — local dissipation equivalence discussed in Section 3. 

5. Results. We present results for the response of the- Kolmogorov flow to an instantaneous change in 
the- amplitude of the forc ing at three different wave numbers. In each case, four runs were made: a direct 
numerical simulation, a large-eddy simulation with the Smagorinskv model with a value for the Smagorinskv 
constant consistent with that predicted by the history-integral model (C,). a large-eddy simulation with 
the plane-averaged history-integral model, a large-eddy simulation with the- dynamic model and. finally, a 
large-eddy simulation with the Smagorinskv model employing a value for the* constant (C rf ) consistent with 
that predicted by the* dynamic model for each case. The direct-numerical simulations were performed with 
a resolution of C4' ! and the large-eddy simulations were performed with a resolution of 32 :l . 

The reason for two separate* Smagorinsky-model simulations at each of the forcing wave* numbers is our 
desire to make the truest possible comparison between the Smagorinskv model and the two more sophisticated 
models considered. Thus, the history-integral-model simulation is most appropriately compared with a 
Smagorinsky-model simulation performed with a Smagorinskv constant equal to that attained by the history- 
integral model m the steady-state limit. Similarly, in view of the circumstance that the Smagorinskv constant 
predicted by the dynamic model varies little in time, it is most appropriate to compare the dynamic-model 
predictions with those of the Smagorinskv model using a value for the Smagorinskv constant fitted as best 
one c an to the time-history <>f the Smagorinskv constant generated by the dynamic model. This ambivalence 
m the Smagorinskv model, with its arbitrary constant, is, of course, one of its more* significant weaknesses, 
and the fact cannot be ignored that it was by using the results of the* dynamic-model calculation that an 
appropriate* value* of the Smagorinsky model was determined without the expensive trial-and-error searc hing 
of parameter spare that, would otherwise be* necessary. 

In attempting to c ompare turbulence simulations of different types with different resolutions, it is nec- 
essary to choose a method of making all simulations start at the same* state. This initial state* should not be* 
too far m advance of the step change in forcing amplitude whose response we want to look at. given the in- 
evitable deviations of the large-eddy simulations from the direct-numerical simulations even for steady-state 
turbulence, and yet, if we take* a 32 :i velocity field from the large-eddy simulations and use it to start the 
G4 f direct, numerical simulations (or vice versa), there is bound to be some period of adjustment while the* 
velocity field orients itself to the new resolution. Both starting the DNS with LES data and starting the 
LES with DNS data was tried in the course of this work and it was found that the* DNS started with LES 
data continued the steady state with almost no discernible* adjustment period, but the* LES started with the 
DNS data experienced an almost 20% drop in turbulent, kinetic energy as soon as the simulation began. As 
a result of these* tests, all comparisons reported in this paper were made by attaining a steady state with 
the* LES code and using the veloc ity field from this run as the initial state for both LES and DNS runs. 

The* data to be* examined in comparing the runs arc* global quantities of the flow: the kinetic energy, the 
dissipation (resolved and subgrid, in the* case of the large-eddy simulations) and the* work done* bv the force* 
cm the* fluid. The three elements of the* overall energy balance are the kinetic energy, A* = 1 /*2 /i >ox (« a + v 2 + 




Fig 5.1. Case I: Energy k versus time. DNS, thin solid line; Smagormsky LES with history-integral-model constant 
dashed line; history-integral LES , dash-dot line; Smagormsky LES with constant generated by dynamic model, thick solid line 
dynamic-model LES, dotted line. 



FlG. r>.2. Case 1: Work done by external force versus time. Line types as in Figure 5.1. 

w 2 )dx , the dissipation, f, and the work done by the force, w f = Re~ l kj J 0 < u > p sin(k/z) dz: 

dk 

(o.l) 

In the case of the large-eddy simulations, both the dissipation of the resolved-scale motion and the dissipation 
of the sub-grid scales (computed in the course of computing the modeled Reynolds stresses) will be examined. 

In Case 1, the forcing wave number is k f = 1 and the initial turbulent state corresponds to Re = 28.2. 
Doubling the amplitude corresponded to changing the Reynolds number to R = 39.9, and it may be seen 
from Figures 5.1 -5.4 that the time histories of the observed global quantities are affected fairly gradually by 
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Fmcj. 5.3. Case 1: Super-grid dissipation versus time . Line types as in Figure 5.1. 



t ic;. 5.1. Case t: Sub-grid dissipation versus time. Line types as in Figure 5.1. 


the abrupt change' in forcing. (The 1 exception is the work done, which jumps simply because it is computed 
directly from the force.) All the LES runs follow the DNS results fairly closely until t ~ 6 or 7, after which 
the runs with the dynamic model and with the Smagorinsky model with the constant suggested by the 
dynamic-model run continue' to represent the DNS results tolerably well, but the history-integral model and 
the Smagorinsky run with C h deviate' significantly from the DNS results. 

In Case 1 2, the forcing wave number is kf =4 and the initial and final Reynolds numbers are' Re = 28.0 
and Re — 39. G. In this case, the response of the* measured global quantities (Figures 5.5 5.8) is more abrupt: 
there* is a rapid rise* in most of the* plotter! quantities when the* forcing is changed. Wo begin to se*e' more 
se'rious discrepane*ies between the LES and DNS runs, and the discrepancies appear earlier. For example. 
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Fig. 5.5. Case 2: Energy k versus time . Line types as in Figure 5.1. 



Fig. 5.6. Case 2: Work done by external force versus time. Line types as in Figure 5.1. 


the history-integral-model run and the associated Smagorinsky-model run fail to attain the proper level of 
kinetic energy very early in the run. The dynamic-model run and its associated Smagorinskv run do correctly 
predict the initial peak in the kinetic energy, but begin to deviate from the DNS results fairly soon afterward. 
The other plotted quantities show similar deviations. 

Case 3 is for a wave number kj — G. The initial and final Reynolds numbers are Re = 28.1 and Re = 39.7. 
Figures 5.9-5.12 show that the deviations noted in Case 2 become even more pronounced in this case: the 
history-integral-model run and its associated Smagorinsky run underestimate the kinetic energy throughout 
the run; the dynamic-model run and its associated Smagorinsky-model run overestimate the kinetic energy 
for most of the run, overshooting the initial peak significantly, as well. 
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These results indicate strongly that the plane-averaged history-integral model provides no significant 
improvement over the Smagorinsky model for this problem. The dynamic model does predict the time his- 
tories of global quantities studied here significantly better than the history-integral model and its associated 
Smagorinsky model, but significant deviations from the direct-numerical simulations are still present, par- 
ticularly for the higher- wavenumber forcing. However, the Smagorinsky-model run whose constant was fixed 
at. the dynamic-model level did nearly as well; this would seem to indicate that the strength of the dynamic 
model for this problem is its ability to automatically choose a good value for the Smagorinsky constant, 
latliei than its ability to accomodate spatial and temporal variations in that ''constant.." 
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Fig. 5.9. Case 3: Eneiyy k versus time. Line types as in Figure 5.1. 



KiG. 5.10. Case 3: Work done by external force versus time. Line types as in Figure 5.1. 

6. Conclusion. The purpose of this final section is to make some progress towards understanding 
the results of the previous section and their consequences for the large-eddy simulation of non-equilibrium 
turbulent flows. Contrasts with the analogous investigation of Nwafor and Woodruff [6] for accelerated plane 
channel flow will also be discussed. 

We are concerned with two basic questions. The more fundamental question is the extent to which, and 
under what conditions, a steady-state model like the Smagorinsky model breaks down for non-equilibrium 
turbulent flows. The secondary question is whether or not the time-dependent history-integral model pro- 
posed in [6] represents any improvement over the Smagorinsky model. 

Considering these questions in turn, it is clear from the results of the previous section that the agreement 
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FlG. 5.1 1 . Case 3: Super-grid dissipation versus time. Line types as in Figure 5.1. 



Fig. 5.12. Case 3: Sub-grid dissipation versus tune. Line types as in Figure 5.1. 


between LES and DNS deteriorates as the wavenumber of the forcing increases. It also seems fair to assort 
that, even in the ease of the best agreement between LES and DNS, that of Case 1 , with kf = 1 , quantities 
based on the mean velocities, such as the bulk velocities, the work and the overall kinetic energy, are 
simulated more accurately by the LES than are quantities determined by the fluctuating velocity field, such 
as the dissipation. The present results thus suggest, that the steady-state Smagorinsky model becomes less 
satisfactory for simulating 11011-equilibrium turbulent flows as the gradients in the flow become larger (a 
criterion we shall attempt to make' more precise below) and as one's interests are concentrated on higher- 
order moments of the flow. 

It is also clear that the plane-averaged history-integral model exhibits no improvement over the Smagorin- 
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Table 6.1 
Summary of Cases. 


Case 

k f 

Initial Re 

Final Re 

INI 

Re f 

1 

1 

28.2 

39.9 

1.00 

2.2 

2 

4 

28.0 

39.6 

0.85 

15 

3 

6 

28.1 

39.7 

0.85 

34 

Ref. [12] 

* 

* 

* 

* 

43 


sky model in the present investigation. In those cases where the results of the LES with the two models do 
deviate, the deviation is on the order of the noise in the signal and is often in the wrong direction, anyway. 
We consequently conclude that the discrepancies between the DNS and the Smagorinsky LES observed in 
the present calculations are not the result of history effects, at least insofar as history effects are represented 
by the plane-averaged history-integral model. 

This second conclusion is in marked contrast to that indicated by the results of the accelerated plane 
channel flow simulations of Nwafor and Woodruff [6]. There, too, LES with the Smagorinsky and the 
plane-averaged history-integral models were compared with DNS for an impulsively-disturbed flow. Time- 
histories of global quantities contributing to the energy balance of the flow were examined (production and 
dissipation), as well as the time history of the wall shear stress. The simulations showed that the history- 
integral model significantly improved the LES predictions for the production and the dissipation in the initial 
stage of the response to the disturbance and qualitatively improved some aspects of the wall shear stress 
predictions. 

In light of all this, it would be desirable to characterize as precisely as possible the differences between 
the three cases of the present investigation which lead to the increasing discrepancy between the LES and 
the DNS, as well as the differences between the present fiow T and the plane channel flow of [6] which cause 
the simulations of one to be improved by the history-integral model but not those of the other. 

Other than the increasing discrepancy between the LES and DNS as the wave number of the forcing 
increases, the most obvious difference between the results of the three cases is the relative magnitude of the 
resolved and subgrid portions of the dissipation. These two quantities may be regarded as characteristic of 
the magnitudes of the viscous and Reynolds-stress terms in the filtered-averaged equations and so their ratio 
is likely to be an important non-dimensional parameter describing the interplay between the resolved and 
sub-grid dynamics. The ratio of local values of the sub- grid to the resolved dissipation is 


(C 8 A)- (SnnSnnY'-SijSjj = cjA (AS mtl AS mn ) 

vSjjSjj v 


A characteristic value of this ratio may be formed with a characteristic strain rate So , 


( 6 . 2 ) 


, A ■ ASo 

r~ 


It may clearly be viewed as a Reynolds number based on the grid size A and the velocity ASo- This last 
velocity scale may be regarded as a characteristic velocity difference induced by the shear for neighboring 
points on the grid. We will consequently define Rc ( = A 2 So/v and examine the success of this parameter in 
predicting the behavior of the simulations. (The Smagorinsky constant c s has been dropped; it is irrelevant 
to our purpose of comparing the relative magnitudes between cases.) 

Let us assume that the magnitude of the shear in the flows studied here may be characterized by the 
magnitude of the shear of the (plane- averaged) mean- velocity profile. Then we may take So = Ay||u||, 
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where ||u|| is the mean peak value of the ^-component of velocity in the initial turbulent state, and so find 
that in the three cases of the present investigation Re< has, just before the instantaneous change in forcing 
amplitude, the values 2.2, 15 and 34, respectively. (See Table G.l) Correlating these values with the results 
described in the previous section, we see that smaller (0(1)) values of Re f seem to correspond to fairly good 
agreement between the LES and DNS, no measurable difference between the LES’s with the history-integral 
model and with the Sinagorinsky model with C/, and much more dissipation in the resolved scales than in 
the sub-grid scales. Intermediate values of Rc t on the order of a dozen or so seem to correspond to larger 
discrepancies between LES and DNS and the resolved dissipation exceeds the sub-grid dissipation by an order 
of magnitude. Finally, values of Re ( on the order of 30 or 40 seem to correspond to large deviations between 
LES and DNS and to a sub-grid dissipation that is commensurate with the dissipation of the resolved scales. 

An additional data point comes from the channel-flow computations of Nwafor and Woodruff [6]: using 
the shear at the wall as the characteristic value of the shear, one finds that Rr ( at tin* initial state has the 
value 43. The disci epancv between the DNS and the Sinagorinsky LES, as well as the differences between the 
Smagoi insky and plane-averaged history-integral model LES s, are consistent with the conclusions drawn in 
the pievious paiagiaph for simulations with moderate values of R (\ . The difference between these ehannel- 
fl° w results and the results of the present study, of course, is that the new model improved the LES results 
in the channel-flow simulations. No such improvement oceured in the present investigation. 

The Reynolds number Re ( thus seems to provide some indication of whether a large-eddy simulation 
using a steady-state model like the Sinagorinsky model will faithfully reproduce at least the global quantities 
in a non-equilibrium turbulent flow. It does not, however, help to explain when the history-integral model 
employed here will offer an improvement for an LES calculation. 
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